*********************************************************
*********************************************************
*                  RELIGIOSITY PROJECT                  *	
*              Descriptives and Regressions			    *
***      		Last modified in Oct 2023        	 ****
*********************************************************
*********************************************************

***********
**# Index
***********

	* Main Tables: 
		* Table 1: Comparing members of Pentecostal churches and Traditional churches
		* Table 2: Human capital, the choice of religious denomination and religiosity: Descriptive evidence
		* Table 3: The effect of the deworming treatment on the choice of religious denomination and religiosity
	
	
	* Appendix Tables:
		* Table A2: Religious practices and beliefs
		* Table A3: Religious practices and beliefs – more detailed classification of churches	
		* Table A4: List of religious denomination reported by the respondents and their classification into Traditional Christian, Pentecostal and Other categories
		* Table A5: Correlations between measures of human capital and the choice of religious denomination - more detailed classification of churches
		* Table A6: Correlations between measures of human capital and components of the Religiosity index
		* Table A7: The effect of the deworming treatment on the individual components of the Religiosity index
		* Table A8: The effect of the deworming treatment on the choice of religious denomination using multinomial logit regression – more detailed classification of churches
		* Table A9: The effect of the deworming treatment on the choice of religious denomination and religiosity, across the four KLPS rounds
		* Table A10: The effect of the deworming treatment on the choice of religious denomination and religiosity among the sub-sample of older respondents, across the four KLPS rounds
		* Table A11: The effect of the deworming treatment on the choice of religious denomination and religiosity among the sub-sample of younger respondents, across the four KLPS rounds
		* Table A12: The effect of the deworming treatment on the choice of religious denomination and religiosity by gender
		* Table A15: Treatment effect on confidence in modern science and medicine
		* Table A16: School Intraclaster Correlation
		* Table A17: Effect of Peer Pastors on the Conversion of Non-Pastor School Mates, KLPS 4 Sample Only
		
		
	* Main Figures
		* Figure 1: Development of the share of respondents affiliated with Traditional Christian, Pentecostal and Other religious denominations
		* Figure 2: The effects of the deworming treatment on the likelihood of being a member of Pentecostal churches, across the four survey rounds
		
		
	* Appendix Figures
		* Figure A1: Development of the share of respondents affiliated with different religious denominations -more detailed classification
		* Figure A2: Development of importance of religion and church attendance
		* Figure A3: The effects of the deworming treatment on the likelihood of being a member of Pentecostal churches (conservative classification), across the four survey rounds
	
		

***********
**# Setup
***********
	clear all
	set more off
	set maxvar 32000

	
	
*********************
**# Define Programs
*********************
	do "$do/Programs.do"
	
	global controls "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced i.interview_year_KLPS"



*****************
**# Main Tables
*****************

**## Table 1: Comparing members of Pentecostal churches and Traditional churches
	use "$data/Religiosity_Crosssection.dta", clear
	keep if inlist(church_class_minimal_KLPS4, 1, 2)
	
	local col_width 200

	preserve
		label define pent_or_trad 	1 "Pentecostal" 2 "Traditional"
		label value church_class_minimal_KLPS4 pent_or_trad
		
		* Panel A - Views about gender roles
		iebaltab s8_30femalemechanic s8_31maledecisions s8_32husbandchores s8_46femalepolitic s8_47femaletradrole s8_48beatwife, ///
			groupvar(church_class_minimal_KLPS4) nonote rowvarlabel  texcolwidth("`col_width' pt") savetex("$tables/panel_a_gender_pentvtrad.tex") replace 
			
		* Panel B - Consumption
		iebaltab s24_1cigarettes s24_2alcohol , groupvar(church_class_minimal_KLPS4) nonote  rowvarlabel texcolwidth("`col_width' pt") ///
			savetex("$tables/panel_b_consumption_pentvtrad.tex") replace

		* Panel C - Political views
		iebaltab s8_15apart12_rallies s8_15bpart12_demonstr s8_15cpart12_discuss s8_17avote s8_24democracy s8_25politics s8_28leaders, ///
			groupvar(church_class_minimal_KLPS4) nonote rowvarlabel  texcolwidth("`col_width' pt") savetex("$tables/panel_c_poli_pentvtrad.tex")  replace
			
		* Panel D - Trust
		iebaltab s8_4trust s8_5trust_own_tribe s8_6trust_oth_tribe s8_7trust_church s8_8trust_oth_church , ///
			groupvar(church_class_minimal_KLPS4) nonote rowvarlabel  texcolwidth("`col_width' pt") savetex("$tables/panel_d_trust_pentvtrad.tex") replace 
	restore
	
	
	
**## Table 2: Human capital, the choice of religious denomination and religiosity: Descriptive evidence
	use "$data/Religiosity_Panel.dta", clear
	
	preserve
		keep if inlist(survey_round,3,4)
		global outcomes_corr	"curr_trad_KLPS curr_ren_KLPS curr_pente_h_KLPS o12_index_KLPS"
		global covariates_corr  "cov_eduattainment_KLPS3 cov_ravens_KLPS3 cov_fedu_KLPS4 cov_medu_KLPS4 cov_earnings_KLPS" 
		global het1				"older"
		global het2				"younger"
		global panelA			"all participants"
		global panelB			"older participants"
		global panelC			"younger participants"
		global tablename		"T2"
		global header			"Tone"
		correlations
	restore
	
	
	
**## Table 3: The effect of the deworming treatment on the choice of religious denomination and religiosity

	* Column 1
	preserve 
		keep if inlist(survey_round,3)
		global outcomes 	    "cov_eduattainment_KLPS3"
		global het1				"older"
		global het2				"younger"
		global panelA			"all participants"
		global panelB			"older participants"
		global panelC			"younger participants"
		global tablename		"T3_edu"
		global header			"Ttwo"
		regressions
	restore
	
	* Column 2
	preserve 
		keep if inlist(survey_round, 2, 3, 4)
		global outcomes 	    "cov_earnings_KLPS"
		global het1				"older"
		global het2				"younger"
		global panelA			"all participants"
		global panelB			"older participants"
		global panelC			"younger participants"
		global tablename		"T3_earn"
		global header			"Ttwo"
		regressions
	restore
	
	
	* Column 3-8
	preserve
		global outcomes 	    "curr_trad_KLPS curr_ren_KLPS curr_pente_h_KLPS renewal_totyr_share_KLPS  renewal_totyr_share_B_KLPS o12_index_KLPS" 
		global het1				"older"
		global het2				"younger"
		global panelA			"all participants"
		global panelB			"older participants"
		global panelC			"younger participants"
		global tablename		"T3"
		global header			"Ttwo"
		regressions
	restore
	

	
*********************
**# Appendix Tables
*********************


**## Table A2: Religious practices and beliefs
	use "$data/Religiosity_Crosssection.dta", clear
	cap erase "$tables/TA2.tex"		
		matrix final = J(3, 11, 0)
		mat rownames final = "Pentecostal" "Traditional Christian" "Other"
		local j=1 
		foreach var of varlist s7_17pastor_KLPS4  s7_19bprophecy_KLPS4 s7_19ahealing_KLPS4 s7_19cdevil_KLPS4 s7_23witnesstrad_KLPS4 ///
			prayintongues_ever_KLPS4 charismaticprac_ever_KLPS4  s7_24saved_KLPS4 s7_25afterlife_KLPS4 prosperity_agree_KLPS4 health_agree_KLPS4 {
			estpost tabstat `var' if  inlist(church_class_minimal_KLPS4,1,2,3), by(church_class_minimal_KLPS4) statistics(mean) //listwise
			forvalues i=1/3 {
				matrix m1 = e(mean)
				local loc = m1[1,`i']
				matrix final[`i',`j']=`loc'
				matrix list final
			}
			local ++j
		}
		
	esttab matrix(final, fmt(2)) using "$tables/TA2.tex", ///
	replace cell("c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11")  compress mgroups("\vspace{-.25cm}\\ \toprule \vspace{-.2cm}  \\ \mainHeaderTAfour", span) 
			
			
**## Table A3: Religious practices and beliefs – more detailed classification of churches	
	use "$data/Religiosity_Crosssection.dta", clear
	cap erase "$tables/TA3.tex"		
		matrix final = J(5, 11, 0)
		mat rownames final =  "Pentecostal churchesConservative" "Pentecostal churchesBroad" "Catholic" "Protestant Anglican" "Other"
		local j=1
		foreach var of varlist s7_17pastor_KLPS4  s7_19bprophecy_KLPS4 s7_19ahealing_KLPS4 s7_19cdevil_KLPS4 s7_23witnesstrad_KLPS4 ///
			prayintongues_ever_KLPS4 charismaticprac_ever_KLPS4 s7_24saved_KLPS4 s7_25afterlife_KLPS4 prosperity_agree_KLPS4 health_agree_KLPS4 {
			estpost tabstat `var' if  inlist(church_class_minimalless_KLPS4,1,2,3,4,5), by(church_class_minimalless_KLPS4) statistics(mean) //listwise
			forvalues i=1/5 {
				matrix m1 = e(mean)
				local loc = m1[1,`i']
				matrix final[`i',`j']=`loc'
				matrix list final
			}
			local ++j
		}

	esttab matrix(final, fmt(2)) using "$tables/TA3.tex", ///
	replace cell("c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11")  compress mgroups("\vspace{-.25cm}\\ \toprule \vspace{-.2cm}  \\ \mainHeaderTAfour", span) 
		
		
		
**## Table A4: List of religious denomination reported by the respondents and their classification into Traditional Christian, Pentecostal and Other categories
	use "$data/Religiosity_Panel.dta", clear

	* Count number of people in each religion
	collapse (count) pupid (first) church_class_minimalless, by(current_religion_KLPS)
	rename pupid count
	gen percent = count*100/18959
	
	* Sort 
	gen sort_descending_count = count*(-1)
	sort sort_descending_count
	drop sort_descending_count
	
	* Panel B - Traditional
	preserve
		keep if inlist(church_class_minimalless, 3, 4)
		export excel using "$tables/religion_dist_trad.xlsx", firstrow(variables) replace
	restore
	
	* Panel C- Penecostal (large)
	preserve
		keep if inlist(church_class_minimalless, 1)
		export excel using "$tables/religion_dist_pent_h.xlsx", firstrow(variables) replace
	restore
	
	* Panel D - Penecostal (small)
	preserve
		keep if inlist(church_class_minimalless, 2)
		export excel using "$tables/religion_dist_pent_l.xlsx", firstrow(variables) replace
	restore
	
	* Panel E - Penecostal (Other)
	preserve
		keep if inlist(church_class_minimalless, 5)
		export excel using "$tables/religion_dist_other.xlsx", firstrow(variables) replace
	restore
	
	* Panel A - All participants
	collapse (sum) count (sum) percent, by(church_class_minimalless)
	export excel using "$tables/religion_dist.xlsx", firstrow(variables) replace
	
	
	
**## Table A5: Correlations between measures of human capital and the choice of religious denomination - more detailed classification of churches
	use "$data/Religiosity_Panel.dta", clear
	preserve	
		keep if inlist(survey_round,3,4)
		global outcomes_corr	"curr_cath_KLPS curr_angl_KLPS curr_pente_h_KLPS curr_pente_l_KLPS"
		global covariates_corr  "cov_eduattainment_KLPS3 cov_ravens_KLPS3 cov_fedu_KLPS4 cov_medu_KLPS4 cov_earnings_KLPS" 
		global het1				""
		global het2				""
		global panelA			"all participants"
		global tablename		"TA5"
		global header			"TAfive"
		correlations
	restore
	
		

**## Table A6: Correlations between measures of human capital and components of the Religiosity index
	use "$data/Religiosity_Panel.dta", clear
	preserve
		keep if inlist(survey_round,3,4)
		global outcomes_corr	"o1a_importance_KLPS o1b_most_imp_identity_KLPS o1c_moreorless_KLPS o2a_attend_regular_KLPS o2b_attend_lastweek_KLPS o2c_donate_hours_KLPS o2d_donate_money_KLPS"
		global covariates_corr  "cov_eduattainment_KLPS3  cov_ravens_KLPS3 cov_fedu_KLPS4 cov_medu_KLPS4 cov_earnings_KLPS"
		global het1				""
		global het2				""
		global panelA			"all participants"
		global tablename		"TA6"
		global header			"TAsix"
		correlations
	restore
		

		
**## Table A7: The effect of the deworming treatment on the individual components of the Religiosity index
	use "$data/Religiosity_Panel.dta", clear

	preserve	
		global outcomes 	   "o1a_importance_KLPS o1b_most_imp_identity_KLPS o1c_moreorless_KLPS o2a_attend_regular_KLPS o2b_attend_lastweek_KLPS o2c_donate_hours_KLPS o2d_donate_money_KLPS"
		global controls  	   "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced i.interview_year_KLPS"
		global het1				"older"
		global het2				"younger"
		global panelA			"all participants"
		global panelB			"older participants"
		global panelC			"younger participants"
		global tablename		"TA7_new"
		global header			"TAseven"
		regressions	 
	restore
	
	
	
**## Table A8: The effect of the deworming treatment on the choice of religious denomination using multinomial logit regression – more detailed classification of churches

	preserve	
		gen   	trad_disentangle 		=.
		replace trad_disentangle=1 if curr_ren_KLPS==1 | curr_oth_d_KLPS==1
		replace trad_disentangle=2 if curr_cath_KLPS==1 
		replace trad_disentangle=3 if curr_angl_KLPS==1 
		
		gen   	renewal_disentangle 	=.
		replace renewal_disentangle=1 if curr_angl_KLPS==1  | curr_cath_KLPS==1  | curr_oth_d_KLPS==1
		replace renewal_disentangle=2 if curr_pente_h_KLPS==1
		replace renewal_disentangle=3 if curr_pente_l_KLPS==1 
		 
		global outcomes 	   "trad_disentangle renewal_disentangle " 
		global het1				"older"
		global het2				"younger"
		global panelA			"all participants"
		global panelB			"older participants"
		global panelC			"younger participants"
		global tablename		"TA8"
		global header			"TAeight"
		
		mlogitreg	 	
		
	restore
		
	
	
**## Table A9: The effect of the deworming treatment on the choice of religious denomination and religiosity across the four KLPS rounds.
	use "$data/Religiosity_Panel.dta", clear
		replace o12_index_KLPS=0 if survey_round==1

		global outcomes 	   "curr_trad_KLPS curr_ren_KLPS curr_pente_h_KLPS  o12_index_KLPS"
		global controls  	   "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced"
		

		forval i=1(1)4 {
			
			cap erase "$tables/Reg_tA9_round`i'.tex"
			estimates clear
		
			foreach outcome of varlist $outcomes {		
				* save percentiles for trimming					
				sum `outcome', d
				
				*control mean & sd - not trimmed
				sum `outcome' if psdp_treat == 0 & survey_round==`i'  [aw=con_weight] //added surveyround
				loc controlMean = r(mean)
				loc controlSD = r(sd)
				
				* Pooled OLS estimator
				eststo: reg `outcome' psdp_treat $controls  [pw=con_weight] if survey_round==`i' ,  cluster(con_psdpsch98) 
				loc treat = _b[psdp_treat]
				estadd scalar control_mean = `controlMean'
				estadd scalar control_SD   = `controlSD'
				estadd loc Controls Yes, replace
				loc teffect = 100*(`treat'/`controlMean')
				estadd scalar t_effect = round(`teffect', .01)	
				
				* saving number of unique pupid observations in regression sample
				preserve
					keep if e(sample)
					bysort pupid (pupid): gen c`outcome' = 1 if _n==1 
					summ c`outcome'
					estadd scalar N_ind = r(N)
				restore	
			}
							
				
			# delimit ;
			esttab using "$tables/Reg_tA9_round`i'.tex", 
			b(%5.3f) se(%5.3f) r2 nogaps star(* 0.1 ** 0.05 *** 0.01)
			nomtitles nocons nolz 
			title(Deworming Treatment and Importance of Religion)
			label se noobs nonumbers collabels(none) replace longtable fragment nomtitles booktabs
			mgroups("\vspace{-.4cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{KLPS `i' - `label'}}} \\ \toprule   \\ \mainHeaderTone", span) 
			keep(psdp_treat) 
			stats(control_mean control_SD t_effect r2 N_ind N, 
			labels("Control Mean" "Control SD" "Treatment Effect (\%)" "R-squared" "Number Individuals" "Number Observations") fmt(2 2 2 2 0))
			substitute("Standard errors in parentheses" " " "\sym{*} \(p<.10\), \sym{**} \(p<.05\), \sym{***} \(p<.01\)" " ")						
			;
			#delimit cr	
			
			estimates clear
		}	

		
	
**## Table A10 and A11: The effect of the deworming treatment on the choice of religious denomination and religiosity across the four KLPS rounds, among sub-samples of younger and older respondents
	use "$data/Religiosity_Panel.dta", clear
		replace o12_index_KLPS=0 if survey_round==1

		global outcomes 	   "curr_trad_KLPS curr_ren_KLPS curr_pente_h_KLPS  o12_index_KLPS"
		global controls  	   "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced"
		
		loc heterogeneity older younger 

		forval i=1(1)4 {							//survey rounds (1 to 4)
			forval j=0(1)1 {							// older vs younger
			
				cap erase "$tables/Reg_tA1011_round`i'_`het'.tex"
				estimates clear
		
				foreach outcome of varlist $outcomes {	
					
					* save percentiles for trimming					
					sum `outcome', d
					
					*control mean & sd - not trimmed
					sum `outcome' if psdp_treat == 0 & survey_round==`i' & older==`j'  [aw=con_weight]
					loc controlMean = r(mean)
					loc controlSD = r(sd)
					
					* Pooled OLS estimator
					eststo: reg `outcome' psdp_treat $controls  [pw=con_weight] if survey_round==`i' & older==`j' ,  cluster(con_psdpsch98) 
					loc treat = _b[psdp_treat]
					estadd scalar control_mean = `controlMean'
					estadd scalar control_SD   = `controlSD'
					estadd loc Controls Yes, replace
					loc teffect = 100*(`treat'/`controlMean')
					estadd scalar t_effect = round(`teffect', .01)	
						
					* saving number of unique pupid observations in regression sample
					preserve
						keep if e(sample)
						bysort pupid (pupid): gen c`outcome' = 1 if _n==1 
						summ c`outcome'
						estadd scalar N_ind = r(N)
					restore	
				}
				
				if `j' == 1 {
					local label "Sample - older participants"
				}
				else if `j' == 0 {
					local label "Sample - younger participants"
				}	
					
					
				# delimit ;
				esttab using "$tables/Reg_tA1011_round`i'_`j'.tex", 
				b(%5.3f) se(%5.3f) r2 nogaps star(* 0.1 ** 0.05 *** 0.01)
				nomtitles nocons nolz 
				title(Deworming Treatment and Importance of Religion)
				label se noobs nonumbers collabels(none) replace longtable fragment nomtitles booktabs
				mgroups("\vspace{-.4cm}\\ \toprule \vspace{-.4cm} \\ &\multicolumn{7}{c}{\emph{\textbf{KLPS `i' - `label'}}} \\ \toprule   \\ \mainHeaderTone", span) 
				keep(psdp_treat) 
				stats(control_mean control_SD t_effect r2 N_ind N, 
						labels("Control Mean" "Control SD" "Treatment Effect (\%)" "R-squared" "Number Individuals" "Number Observations") fmt(2 2 2 2 0))
				substitute("Standard errors in parentheses" " "
										"\sym{*} \(p<.10\), \sym{**} \(p<.05\), \sym{***} \(p<.01\)" " ")						
				;
				#delimit cr	
				estimates clear
			}				
		}
	
	use "$data/Religiosity_Panel.dta", clear 
	
	
	
**## Table A12: The effect of the deworming treatment on the choice of religious denomination and religiosity by gender
	preserve
		* Gen male dummy
		gen male = female + 1
		recode male 2=0
		
		global outcomes 	    "curr_trad_KLPS curr_ren_KLPS curr_pente_h_KLPS renewal_totyr_share_KLPS  renewal_totyr_share_B_KLPS o12_index_KLPS" 
		global controls  	    "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced i.interview_year_KLPS"
		global het1				"female"
		global het2				"male"
		global panelA			"all participants"
		global panelB			"female sample"
		global panelC			"male sample"
		global tablename		"TA12"
		global header			"Ttwo"
		regressions
	restore
	
	


	
	
	
**## Table A15: Treatment effect on confidence in modern science and medicine
	use "$data/Religiosity_Panel.dta", clear
	keep if survey_round == 4
	
	preserve	
		global outcomes 	   "child_bednet_KLPS4 child_vaccine_KLPS4"
		global controls  	   "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced i.interview_year_KLPS"
		global het1				""
		global het2				""
		global tablename		"TA15"
		global header			"TAfifteen"
		regressions	 
	restore
	
	
	
**## Table A16: School Intraclaster Correlation
	use "$data/Religiosity_Crosssection.dta", clear
	
	* Ever converted from Trad. to Pent.
	loneway trad_to_ren con_psdpsch98
	
	* Ever converted to Pent.
	loneway converted_to_ren con_psdpsch98
	
	* Currently Pent. 
	loneway curr_ren con_psdpsch98	
	
	
**## Table A17: Effect of Peer Pastors on the Conversion of Non-Pastor School Mates, KLPS 4 Sample Only
	use "$data/Religiosity_Crosssection.dta", clear
	
	preserve
		keep con_psdpsch98 s7_17pastor
		collapse (count) s7_17pastor , by(con_psdpsch98)
		rename s7_17pastor num_pastor_klps4
		save "$data/KLPS4_pastor.dta", replace
	restore
	
	merge m:1 con_psdpsch98 using "$data/KLPS4_pastor.dta", nogen keepusing(num_pastor_klps4)
	erase "$data/KLPS4_pastor.dta"
	
	
	gen trad_to_ren_KLPS4 = 0
	gen converted_to_ren_KLPS4 = 0
	
	foreach i of numlist 2015/2020{
		local j = `i' + 1
		replace trad_to_ren_KLPS4 = 1 if renewal_min_`i' == 0 & renewal_min_`j' == 1 & trad_min_`i' == 1 & trad_min_`j' == 0
		replace converted_to_ren_KLPS4 = 1 if renewal_min_`i' == 0 & renewal_min_`j' == 1
	}
	label var trad_to_ren_KLPS4 "Ever converted from Trad. to Pent."
	label var converted_to_ren_KLPS4 "Ever converted to Pent."

	
	label var num_pastor_klps4 "Number of Pastor Cohorts"

	global controls "con_cost_sharing con_saturation_dm_KLPS con_demeaned_popT_6k_KLPS con_zoneidI2-con_zoneidI8 con_pup_pop_KLPS con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2-con_std98_base_I6 con_avgtest96 con_wave2 female voced"
	

	eststo clear
	eststo: quietly reg trad_to_ren_KLPS4 			num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0, vce(cluster con_psdpsch98)
	eststo: quietly reg converted_to_ren_KLPS4 		num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0, vce(cluster con_psdpsch98)
	eststo: quietly reg curr_ren_KLPS4 				num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0, vce(cluster con_psdpsch98)
	eststo: quietly reg curr_trad_KLPS4 			num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0, vce(cluster con_psdpsch98)
	esttab using "$tables/pastor_overall.tex", label se replace mtitle("Ever converted from Trad. to Pent." "Ever converted to Pent." "Pent. Christain" "Trad. Christain") ///
	nobaselevels booktabs nolz style(tex) compress keep(num_pastor_klps4) nonotes
	eststo clear
	
	eststo clear
	eststo: quietly reg trad_to_ren_KLPS4 			num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor==0 & older==1, vce(cluster con_psdpsch98)
	eststo: quietly reg converted_to_ren_KLPS4 		num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor==0 & older==1, vce(cluster con_psdpsch98)
	eststo: quietly reg curr_ren_KLPS4 				num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0 & older==1, vce(cluster con_psdpsch98)
	eststo: quietly reg curr_trad_KLPS4 			num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0 & older==1, vce(cluster con_psdpsch98)
	esttab using "$tables/pastor_older.tex", label se replace mtitle("Ever converted from Trad. to Pent." "Ever converted to Pent." "Pent. Christain" "Trad. Christain") ///
	nobaselevels booktabs nolz style(tex) compress keep(num_pastor_klps4) nonotes nonumbers
	eststo clear
	
	eststo clear
	eststo: quietly reg trad_to_ren_KLPS4 			num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor==0 & older==0, vce(cluster con_psdpsch98)
	eststo: quietly reg converted_to_ren_KLPS4 		num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor==0 & older==0, vce(cluster con_psdpsch98)
	eststo: quietly reg curr_ren_KLPS4 				num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0 & older==0, vce(cluster con_psdpsch98)
	eststo: quietly reg curr_trad_KLPS4 			num_pastor_klps4 $controls [pw=con_weight_KLPS] if s7_17pastor ==0 & older==0, vce(cluster con_psdpsch98)
	esttab using "$tables/pastor_younger.tex", label se replace mtitle("Ever converted from Trad. to Pent." "Ever converted to Pent." "Pent. Christain" "Trad. Christain") ///
	nobaselevels booktabs nolz style(tex) compress keep(num_pastor_klps4) nonumbers
	eststo clear	
	
	
	
	
	
	
	
******************
**# Main Figures
******************
use "$data/Religiosity_Crosssection.dta", clear
	
	
**## Figure 1: Development of the share of respondents affiliated with Traditional Christian, Pentecostal and Other religious denominations
	preserve
		greshape long 	renewal_min_ trad_min_ oth_min_ , i(pupid) j(year)
		gcollapse (mean) renewal_min_  trad_min_  oth_min_ , by(year)	
		
		twoway (connected   renewal_min_ year,  ytitle("", height(5))  xtitle(" ", height(6))  mcolor(midblue) 	msize(medsmall) msymbol(S)	lcolor(midblue) ///
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue) 	lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank))) ///
			   (connected   trad_min_ year,  	ytitle("", height(5))  xtitle(" ", height(6))  mcolor(cranberry) 	msize(medsmall) msymbol(O)	lcolor(cranberry) ///
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue)  lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank))) ///
			   (connected   oth_min_ year,  	ytitle("", height(5))  xtitle(" ", height(6))  mcolor(green)  		msize(medlarge) msymbol(X)	lcolor(green) ///
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue) 	lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(vvvthick) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank)) ///
			   text(.20 1998.9  "22.3%", place(w) size(vsmall)) /// 
			   text(.69 1998.9  "72.2%", place(w) size(vsmall)) /// 
			   text(.03 1998.7  "5.5%", place(w) size(vsmall)) /// 
			   text(.53 2020.9  "50.1%", place(w) size(vsmall)) /// 
			   text(.40 2020.9  "42.4%", place(w) size(vsmall)) /// 
			   text(.05 2021  "7.5%", place(w) size(vsmall)) /// 
			   text(.08 2003.8  "Other religious denominations", place(w) size(vsmall)) /// 
			   text(.75 2003.5  "Traditional Christian Churches", place(w) size(vsmall)) /// 
			   text(.3 2001.8  "Pentecostal Churches" "(broad classification)", place(w) size(vsmall))), /// 
			   xlabel(, labsize(small) nogrid)  ylabel(, labsize(small)) scheme(s2mono)	legend(off)
		
		  graph export "$tables/f1_yr.pdf", replace 
	restore
	

	
	
**## Figure 2: The effects of the deworming treatment on the likelihood of being a member of Pentecostal churches, across the four survey rounds
	use "$data/Religiosity_Panel.dta", clear

	label var   curr_ren_KLPS 				"Pentecostal (broad classification)"
	rename 		curr_ren_KLPS 				curr_pent
	rename 		older 						o
	rename 		younger 					y

	global outcomes 	    "curr_pent" 
	global het1				"o"
	global het2				"y"
	
	estimates clear
	
	* Calculate statistics 
	foreach outcome of varlist $outcomes {
		global controls "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced i.interview_year_KLPS"
		
		* All FR, Round Pooled
	
			*control mean & sd - not trimmed
			sum `outcome' if psdp_treat == 0 [aw=con_weight] 
			loc controlMean = r(mean)
			loc controlSD = r(sd)
			
			* Pooled OLS estimator
			capture eststo: capture reg `outcome' psdp_treat $controls  [pw=con_weight] ,  cluster(con_psdpsch98) 
			reg `outcome' psdp_treat $controls  [pw=con_weight] ,  cluster(con_psdpsch98) 
			estimate store pooled_all_`outcome'
			loc treat = _b[psdp_treat]
		
			loc controlmean_pooled_all 	= `controlMean'
			loc treat_pooled_all 		= `treat'  
			loc teffect_pooled_all 		= 100*(`treat_pooled_all'/`controlmean_pooled_all')
			lincomest psdp_treat
			eststo pooled_all	
			matrix list r(table)
			matrix define pooled_all_`outcome' = r(table)
	
		* All FR, By Round
		forval i=1(1)4 {
		
			global controls_round "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced"
			
			if "`outcome'"=="o12_index_KLPS" & survey_round==1 {
				sum survey_round
			}
			else {
				*control mean & sd - not trimmed
				sum `outcome' if psdp_treat == 0 & survey_round==`i' [aw=con_weight] 
				loc controlMean = r(mean)
				loc controlSD = r(sd)
				
				* Pooled OLS estimator
				capture eststo:  reg `outcome' psdp_treat $controls_round if survey_round==`i' [pw=con_weight] ,  cluster(con_psdpsch98) 
				reg `outcome' psdp_treat $controls_round if survey_round==`i' [pw=con_weight] ,  cluster(con_psdpsch98) 
				capture estimate store klps`i'_`outcome'
				capture loc treat = _b[psdp_treat]	
				
				capture loc controlmean_pooled_all 	= `controlMean'
				capture loc treat_pooled_all 		= `treat'  
				capture loc teffect_pooled_all 		= 100*(`treat_pooled_all'/`controlmean_pooled_all')
				capture lincomest psdp_treat
				capture eststo pooled_all	
				capture matrix list r(table)
				capture matrix define klps`i'_`outcome' = r(table)		
			}		
		}
		
		* By FR Age at Baseline
		loc heterogeneity $het1 $het2 
		foreach het of local heterogeneity   {
			
			* Round Pooled
			
				*control mean & sd - not trimmed
				sum `outcome' if psdp_treat == 0 & `het' == 1 [aw=con_weight] 
				loc controlMean = r(mean)
				loc controlSD = r(sd)
				
				* Pooled OLS estimator by gender
				eststo: reg `outcome' psdp_treat $controls  if `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
				reg `outcome' psdp_treat $controls  if `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
				estimate store pooled_all_`outcome'_`het'
				loc treat = _b[psdp_treat]
				
			* By Round
			forval i=1(1)4 {
				if "`outcome'"=="index" & survey_round==1 {
					sum survey_round
				}
				else {
					*control mean & sd - not trimmed
						sum `outcome' if psdp_treat == 0 & survey_round==`i' [aw=con_weight] 
						loc controlMean = r(mean)
						loc controlSD = r(sd)
					* Pooled OLS estimator
						capture eststo:  reg `outcome' psdp_treat $controls_round if survey_round==`i' & `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
						reg `outcome' psdp_treat $controls_round if survey_round==`i' & `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
						capture estimate store klps`i'_`outcome'_`het'
				}		
			}
		}
	}

	* Draw Figure
	foreach outcome of varlist $outcomes {	

		local lbl: var label `outcome'
		coefplot (pooled_all_`outcome', msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(0)) ///
				(pooled_all_`outcome'_o, 	msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(0.5)) ///
				(pooled_all_`outcome'_y,msymbol(T)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(1)) 	///
				(klps1_`outcome', 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(3.5)) ///
				(klps1_`outcome'_o, 		msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(4)) ///
				(klps1_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(4.5)) 	///	
				(klps2_`outcome',	 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(6)) ///
				(klps2_`outcome'_o,			msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(6.5)) ///
				(klps2_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(7)) 	///	
				(klps3_`outcome', 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(8.5)) ///
				(klps3_`outcome'_o, 		msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(9)) ///
				(klps3_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(9.5)) 	///	
				(klps4_`outcome', 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(11)) ///
				(klps4_`outcome'_o, 		msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(11.5)) ///
				(klps4_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(12)) 	///
			, yline(0, lcolor(black)) xline(3.3, lcolor(bluishgray ) lpatter(dash))	 xlabel(1.5 `" "Pooled" "(1998-2021)" "' 5 `" "KLPS1" "(2003-2005)" "' ///
			7.5 `" "KLPS2" "(2007-2009)" "' 10 `" "KLPS3" "(2011-2014)" "' 12.5 `" "KLPS4" "(2017-2021)" "' 13.5 " ", labsize(small) noticks)  /// 
			ylabel(-.15  -.1 -.05 0 .05 0.1 , labsize(small)) relocate(pooled_all_`outcome'=1.6)  vertical keep(psdp_treat)	///
			plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank)) ///
			graphregion(fcolor(white) lcolor(none) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank)) ///
			ytitle("Deworming treatment effects", size(small)) ysc(titlegap(3)) title(" ", position(6) margin(b=8)) ///
			legend(order( 2 "All" 4 "Older than 12" 6 "Younger than 12") col(8) nobox region(style(none))) legend( ring(0) position(12))  ///
			text( -.22 0 	"Control" "Mean", size(small))  ///
			text( -.22 1.5 	"[All:0.40]" "[Old:0.44]" "[Young:0.37]", size(small))  /// *pooled
			text( -.22 5 	"[All:0.29]" "[Old:0.35]" "[Young:0.26]", size(small))  /// *KLPS1
			text( -.22 7.5 	"[All:0.36]" "[Old:0.40]" "[Young:0.34]", size(small))  /// *KLPS2
			text( -.22 10 	"[All:0.45]" "[Old:0.50]" "[Young:0.41]", size(small))  /// *KLPS3
			text( -.22 12.5 "[All:0.52]" "[Old:0.55]" "[Young:0.50]", size(small))   //*KLPS4

		graph export "$tables/plot_`outcome'.pdf", replace 
	}
	



**********************
**# Appendix Figures
**********************	

**## Figure A1: Development of the share of respondents affiliated with different religious denominations - More detailed classification
	use "$data/Religiosity_Crosssection.dta", clear
	preserve
		greshape long 	pente_l_min_ pente_h_min_ cath_min_ angl_min_ oth_min_ , i(pupid) j(year)
		gcollapse (mean) pente_l_min_ pente_h_min_ cath_min_ angl_min_ oth_min_ , by(year)	

		twoway (connected   pente_l_min_ year, ytitle("", height(5)) xtitle(" ", height(6))  mcolor(cranberry) 		msize(medsmall) msymbol(S) lcolor(cranberry) /// 	
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue) lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank))) ///
			   (connected   pente_h_min_ year, ytitle("", height(5)) xtitle(" ", height(6))  mcolor(pink) 			msize(medsmall) msymbol(O) lcolor(pink) ///
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue) lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank))) ///
			   (connected   cath_min_ year, 	   ytitle("", height(5)) xtitle(" ", height(6)) mcolor(midblue) 		msize(medsmall) msymbol(D) lcolor(midblue) 	///
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue) lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank))) ///
			   (connected   angl_min_ year,	   ytitle("", height(5)) xtitle(" ", height(6)) mcolor(eltblue) 		msize(medlarge) msymbol(X) lcolor(eltblue) 	///
					lwidth(medthick) graphregion(fcolor(white) lcolor(midblue) lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank))) ///
			   (connected   oth_min_ year,	   ytitle("", height(5)) xtitle(" ", height(6)) mcolor(green) 			msize(medlarge) msymbol(|) lcolor(green) ///
					graphregion(fcolor(white) lcolor(midblue) lwidth(none) lpattern(blank) ifcolor(none) ilcolor(none) ilwidth(none) ilpattern(blank)) ///
					plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank)) ///
			   text(.196 	2000.99  "Large Pentecostal" "Churches", place(w) size(vsmall)) /// 
			   text(.104	2000.99  "Small Pentecostal" "Churches", place(w) size(vsmall)) /// 
			   text(.473 	1999.2  "Catholic", place(w) size(vsmall)) /// 
			   text(.285 	2001.2  "Protestant Anglican", place(w) size(vsmall)) /// 
			   text(.04 	2003.2  "Other religious denominations", place(w) size(vsmall))), /// 
			   xlabel(, labsize(small) nogrid)  ylabel(, labsize(small)) scheme(s2mono)	legend(off)

		graph export "$tables/fa1_yr.pdf", replace
	restore	
	
	
	
**## Figure A2: Development of importance of religion and church attendance
	use "$data/Religiosity_Panel.dta", clear
	preserve
	
		forval i=1(1)4 {
			gen o1a_importance_KLPS`i'=o1a_importance_KLPS if survey_round==`i'
			g 		relvimp`i'=1 if o1a_importance_KLPS`i'==3
			replace relvimp`i'=0 if inlist(o1a_importance_KLPS`i',2,1)
			
			gen o2b_attend_lastweek_KLPS`i'= o2b_attend_lastweek_KLPS if survey_round==`i'
		}
		
		collapse (mean) relvimp* o2b_attend_lastweek_KLPS*
		
		g id=1
		drop o2b_attend_lastweek_KLPS
		
		reshape long relvimp o2b_attend_lastweek_KLPS, j(round) i(id)
		
		label var relvimp 					"Self-reported importance (very important)"
		label var o2b_attend_lastweek_KLPS 	"Attended church last week"
		
		twoway (connected  relvimp round, lcolor(midblue) msymbol(S) mcolor(midblue)) ///
				(connected o2b_attend_lastweek_KLPS round, lcolor(cranberry) msymbol(D) mcolor(cranberry)), ///
				ylabel(0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1) xlabel(0.8 " " 1 `" "KLPS1" "(2003-2005)" "' 2 `" "KLPS2" "(2007-2009)" "' 3 `" "KLPS3" "(2011-2014)" "' ///
				4 `" "KLPS4" "(2017-2021)" "' 4.3 " ", labsize(small) noticks) xtitle("") legend(row(2) position(6)) ///
				graphregion(color(white)) bgcolor(white)
				
		 graph export "$tables/trends_fX.pdf", replace
	
	restore	
	
	
	
**## Figure A3: The effects of the deworming treatment on the likelihood of being a member of Pentecostal churches (conservative classification), across the four survey rounds

	use "$data/Religiosity_Panel.dta", clear


	label var   curr_pente_h_KLPS 			"Pentecostal (conservative classification)"
	rename 		curr_pente_h_KLPS 			curr_pent_h
	rename 		older 						o
	rename 		younger 					y
	global outcomes 	    "curr_pent_h " 
	global het1				"o"
	global het2				"y"
	

	estimates clear
	
	* Calculate statistics 
	foreach outcome of varlist $outcomes {
		global controls "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced i.interview_year_KLPS"
		
		* All FR, Round Pooled
	
			*control mean & sd - not trimmed
			sum `outcome' if psdp_treat == 0 [aw=con_weight] 
			loc controlMean = r(mean)
			loc controlSD = r(sd)
			
			* Pooled OLS estimator
			capture eststo: capture reg `outcome' psdp_treat $controls  [pw=con_weight] ,  cluster(con_psdpsch98) 
			reg `outcome' psdp_treat $controls  [pw=con_weight] ,  cluster(con_psdpsch98) 
			estimate store pooled_all_`outcome'
			loc treat = _b[psdp_treat]
		
			loc controlmean_pooled_all 	= `controlMean'
			loc treat_pooled_all 		= `treat'  
			loc teffect_pooled_all 		= 100*(`treat_pooled_all'/`controlmean_pooled_all')
			lincomest psdp_treat
			eststo pooled_all	
			matrix list r(table)
			matrix define pooled_all_`outcome' = r(table)
	
		* All FR, By Round
		forval i=1(1)4 {
		
			global controls_round "con_cost_sharing con_saturation_dm con_demeaned_popT_6k con_zoneidI2-con_zoneidI8 con_pup_pop con_month_interviewI1_KLPS con_month_interviewI2_KLPS con_month_interviewI3_KLPS con_month_interviewI4_KLPS con_month_interviewI5_KLPS con_month_interviewI6_KLPS con_month_interviewI7_KLPS con_month_interviewI8_KLPS con_month_interviewI9_KLPS con_month_interviewI10_KLPS con_month_interviewI11_KLPS con_month_interviewI12_KLPS  con_std98_base_I2 -con_std98_base_I6 con_avgtest96 con_wave2  female voced"
			
			if "`outcome'"=="o12_index_KLPS" & survey_round==1 {
				sum survey_round
			}
			else {
				*control mean & sd - not trimmed
				sum `outcome' if psdp_treat == 0 & survey_round==`i' [aw=con_weight] 
				loc controlMean = r(mean)
				loc controlSD = r(sd)
				
				* Pooled OLS estimator
				capture eststo:  reg `outcome' psdp_treat $controls_round if survey_round==`i' [pw=con_weight] ,  cluster(con_psdpsch98) 
				capture reg `outcome' psdp_treat $controls_round if survey_round==`i' [pw=con_weight] ,  cluster(con_psdpsch98) 
				capture estimate store klps`i'_`outcome'
				capture loc treat = _b[psdp_treat]	
				
				capture loc controlmean_pooled_all 	= `controlMean'
				capture loc treat_pooled_all 		= `treat'  
				capture loc teffect_pooled_all 		= 100*(`treat_pooled_all'/`controlmean_pooled_all')
				capture lincomest psdp_treat
				capture eststo pooled_all	
				capture matrix list r(table)
				capture matrix define klps`i'_`outcome' = r(table)		
			}		
		}
		
		* By FR Age at Baseline
		loc heterogeneity $het1 $het2 
		foreach het of local heterogeneity   {
			
			* Round Pooled
			
				*control mean & sd - not trimmed
				sum `outcome' if psdp_treat == 0 & `het' == 1 [aw=con_weight] 
				loc controlMean = r(mean)
				loc controlSD = r(sd)
				
				* Pooled OLS estimator by gender
				eststo: reg `outcome' psdp_treat $controls  if `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
				reg `outcome' psdp_treat $controls  if `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
				estimate store pooled_all_`outcome'_`het'
				loc treat = _b[psdp_treat]
				
			* By Round
			forval i=1(1)4 {
				if "`outcome'"=="index" & survey_round==1 {
					sum survey_round
				}
				else {
					*control mean & sd - not trimmed
						sum `outcome' if psdp_treat == 0 & survey_round==`i' [aw=con_weight] 
						loc controlMean = r(mean)
						loc controlSD = r(sd)
					* Pooled OLS estimator
						capture eststo:  reg `outcome' psdp_treat $controls_round if survey_round==`i' & `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
						capture reg `outcome' psdp_treat $controls_round if survey_round==`i' & `het'==1 [pw=con_weight] ,  cluster(con_psdpsch98) 
						capture estimate store klps`i'_`outcome'_`het'
				}		
			}
		}
	}

	* Draw Figure
	foreach outcome of varlist $outcomes {	

		local lbl: var label `outcome'
		coefplot (pooled_all_`outcome', msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(0)) ///
				(pooled_all_`outcome'_o, 	msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(0.5)) ///
				(pooled_all_`outcome'_y,msymbol(T)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(1)) 	///
				(klps1_`outcome', 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(3.5)) ///
				(klps1_`outcome'_o, 		msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(4)) ///
				(klps1_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(4.5)) 	///	
				(klps2_`outcome',	 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(6)) ///
				(klps2_`outcome'_o,			msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(6.5)) ///
				(klps2_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(7)) 	///	
				(klps3_`outcome', 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(8.5)) ///
				(klps3_`outcome'_o, 		msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(9)) ///
				(klps3_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(9.5)) 	///	
				(klps4_`outcome', 	msymbol(D) mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(cranberry) ciopts(color(cranberry)) offset(11)) ///
				(klps4_`outcome'_o, 		msymbol(S)  mlabel(@b) mlabsize(vsmall) mlabcolor(black) mlabformat("%9.3f") mcolor(midblue) ciopts(color(midblue)) offset(11.5)) ///
				(klps4_`outcome'_y, 	msymbol(T)  mlabel(@b) mlabsize(vsmall)  mlabcolor(black) mlabformat("%9.3f") mcolor(green) ciopts(color(green)) offset(12)) 	///
			, yline(0, lcolor(black)) xline(3.3, lcolor(bluishgray ) lpatter(dash))	 xlabel(1.5 `" "Pooled" "(1998-2021)" "' 5 `" "KLPS1" "(2003-2005)" "' ///
			7.5 `" "KLPS2" "(2007-2009)" "' 10 `" "KLPS3" "(2011-2014)" "' 12.5 `" "KLPS4" "(2017-2021)" "' 13.5 " ", labsize(small) noticks)  /// 
			ylabel(-.15  -.1 -.05 0 .05 0.1 , labsize(small)) relocate(pooled_all_`outcome'=1.6)  vertical keep(psdp_treat)	///
			plotregion(fcolor(white) lcolor(white) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank)) ///
			graphregion(fcolor(white) lcolor(none) lwidth(none) lpattern(blank) ifcolor(white) ilcolor(white) ilwidth(none) ilpattern(blank)) ///
			ytitle("Deworming treatment effects", size(small)) ysc(titlegap(3)) title(" ", position(6) margin(b=8)) ///
			legend(order( 2 "All" 4 "Older than 12" 6 "Younger than 12") col(8) nobox region(style(none))) legend( ring(0) position(12)) ///
			text( -.22 0 	"Control" "Mean", size(small))  ///
			text( -.22 1.5 	"[All:0.25]" "[Old:0.29]" "[Young:0.23]", size(small))  /// *pooled
			text( -.22 5 	"[All:0.19]" "[Old:0.22]" "[Young:0.17]", size(small))  /// *KLPS1
			text( -.22 7.5 	"[All:0.20]" "[Old:0.23]" "[Young:0.18]", size(small))  /// *KLPS2
			text( -.22 10 	"[All:0.27]" "[Old:0.31]" "[Young:0.24]", size(small))  /// *KLPS3
			text( -.22 12.5 "[All:0.39]" "[Old:0.42]" "[Young:0.36]", size(small))  // *KLPS4

		graph export "$tables/plot_`outcome'.pdf", replace 
	}	
	
	
	
	
	
**# Statistics in Paper

	use "$data/Religiosity_Panel.dta", clear

	* 95% say religion is very important in their lives. (p11)
	tab o1a_importance_KLPS
	/*
	o1a_importance_KLP |
					 S |      Freq.     Percent        Cum.
	-------------------+-----------------------------------
	Not very important |         95        0.51        0.51
	Somewhat important |        820        4.41        4.92
		Very important |     17,686       95.08      100.00
	-------------------+-----------------------------------
				 Total |     18,601      100.00

	*/
	
	
	
	* 75% report attending church regularly
	tab o2a_attend_regular_KLPS
	/*
	o2a_attend_ |
	regular_KLP |
			  S |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      4,571       24.57       24.57
			Yes |     14,035       75.43      100.00
	------------+-----------------------------------
		  Total |     18,606      100.00

	*/
	
	
	* Respondents spend a non-negligible share of their household budget (4%) on church
	sum o2dbis_donate_hh_money_KLPS
	/*
	    Variable |        Obs        Mean    Std. dev.       Min        Max
	-------------+---------------------------------------------------------
	o~hh_money~S |      3,318    4.611337    49.50941          0    1916.26
	*/
	
	
	
	* Share of "Other" religion remonination is 6% (p17)
	tab curr_oth_d_KLPS
	/*
	curr_oth_d_ |
		   KLPS |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |     17,771       93.74       93.74
			  1 |      1,187        6.26      100.00
	------------+-----------------------------------
		  Total |     18,958      100.00
	*/
	
	
	
	use "$data/Religiosity_Crosssection.dta", clear	
	
	* Conversions within Trad. is 7% (P18)
	tab o4i_convert_within_trad
	/*
	Outcome 4i: |
	  Converted |
		 within |
	Traditional |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      6,292       93.35       93.35
			  1 |        448        6.65      100.00
	------------+-----------------------------------
		  Total |      6,740      100.00
	*/
	
	
	
	* Conversion within Pent. is 25% (P18)
	tab o4i_convert_within_ren
	/*
	Outcome 41: |
	  Converted |
		 within |
	Pentecostal |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      5,080       75.37       75.37
			  1 |      1,660       24.63      100.00
	------------+-----------------------------------
		  Total |      6,740      100.00
	*/
	
	
	use "$data/Religiosity_Panel.dta", clear
	gen curr_trad_KLPS1 = curr_trad_KLPS if survey_round == 1
	gen curr_trad_KLPS4 = curr_trad_KLPS if survey_round == 4
	gen curr_ren_KLPS1 = curr_ren_KLPS if survey_round == 1
	gen curr_ren_KLPS4 = curr_ren_KLPS if survey_round == 4
	
	* Share of Trad. dropped from 65% in KLPS1 to 42% in KLPS4 (P18)
	tab curr_trad_KLPS1
	tab curr_trad_KLPS4
	/*
	curr_trad_K |
		   LPS1 |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      1,802       34.92       34.92
			  1 |      3,359       65.08      100.00
	------------+-----------------------------------
		  Total |      5,161      100.00
		  

	curr_trad_K |
		   LPS4 |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      2,439       57.65       57.65
			  1 |      1,792       42.35      100.00
	------------+-----------------------------------
		  Total |      4,231      100.00
	*/
	
	
	* Share of Trad. dropped from 29% in KLPS1 to 50% in KLPS4 (P18)
	tab curr_ren_KLPS1
	tab curr_ren_KLPS4
	/*
	curr_ren_KL |
			PS1 |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      3,669       71.09       71.09
			  1 |      1,492       28.91      100.00
	------------+-----------------------------------
		  Total |      5,161      100.00

	curr_ren_KL |
			PS4 |      Freq.     Percent        Cum.
	------------+-----------------------------------
			  0 |      2,112       49.92       49.92
			  1 |      2,119       50.08      100.00
	------------+-----------------------------------
		  Total |      4,231      100.00
	*/

	* The share of being "Other" is consistent between 5% to 8% (P18)
	foreach i of numlist 1/4{
		sum curr_oth_d_KLPS if survey_round == `i'
	}

	

**## Percentage of Conversions

	use "$data/Religiosity_Crosssection.dta", clear
	
	* Always trad and always pentecostal
	gen always_trad = 1
	gen always_ren 	= 1
	foreach i of numlist 1998/2020{
		local j = `i' + 1
		replace always_trad = 0 if trad_min_`i' == 0
		replace always_ren = 0 if renewal_min_`i' == 0
	}

	* 38% of FRs converted from Traditional Christian to Pentecostal churches, or from Pentecostal to Traditional Christian churches
	tab1 o4cmin_convert_never_ren
	
	
	* 43% always belonging to a Traditional Christian church and 18% to a Pentecostal denomination
	tab1 always_trad always_ren

	
	* Changes within the group of Traditional Christian denomination 
	* From Catholic to Anglican or vice versa, at 7% of the sample. Among Pentecostal churches 25 %
	tab1 o4i_convert_within_trad o4i_convert_within_ren
	
